Firing Time Statistics for Driven Neuron Models: Analytic Expressions versus 

Numerics 



Michael Schindler, Peter Talkner, and Peter Hanggi 
Institut fur Physik, Universitat Augsburg, Universitatsstrafie 1, D-86135 Augsburg, Germany 

(Dated: February 4, 2008) 

Analytical expressions are put forward to investigate the forced spiking activity of abstract neuron 
models such as the driven leaky integrate-and-fire model. The method is valid in a wide parameter 
regime beyond the restraining limits of weak driving (linear response) and/or weak noise. The novel 
approximation is based on a discrete state Markovian modeling of the full long-time dynamics with 
time- dependent rates. The scheme yields excellent agreement with numerical Langevin and Fokker- 
Planck simulations of the full non-stationary dynamics, not only for the first-passage time statistics, 
but also for the important interspike interval (residence time) distribution. 
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The detailed modeling of neural behavior presents a 
prominent challenge on the intriguing vista towards the 
understanding of neural coding principles. The leaky 
integrate-and-fire (LIF) model, whose deterministic for- 
mulation has been introduced long ago 0, likely is one 
of the most studied abstract neuron models 2]. It is 
characterized by marked simplicity and lack of memory: 
Whenever the neuron has been excited to fire a pulse it 
is reset to a predefined state. The beneficial role of an 
appreciable dose of noise has proved to bestow a key part 
to the interspike statistics of neurons 0]. There exist by 
now numerous studies and important generalizations of 
realistic synaptic models, mostly of a numerical nature, 
that demonstrate the rich behavior of the renewal firing 
probability, e.g. see in Ref. (UOHH- 

In presence of a time-dependent input-stimulation the 
stochastic firing process becomes non-stationary, which 
in turn significantly complicates the stochastic firing 
statistics. Nevertheless, the signal transmission and its 
detection can exhibit a remarkable improvement via the 
phenomenon of Stochastic Resonance jfj. The dynamics 
of the neuronal firing probability emerges due to a large 
bombardment of synaptic spike events; consequently, it 
is customary to employ a diffusion approximation for the 
stochastic dynamics of the membrane potential x(t). The 
complexity of the driven abstract LIF model thus as- 
sumes the archetype, non-stationary Langevin dynamics 

x(t)=-Xx{t)+n + f(t)+V2D€(t) (1) 

where the process starts at a time s at x(s) = xq 
and fires when it reaches the threshold voltage x = a. 
Here, f(t) presents a general, time-dependent stimu- 
lus which, for example, can be chosen to be oscilla- 
tory, and £(t) is white Gaussian noise. The dynamics 
of the process x(t) is equivalently described by a Fokker- 
Planck (FP) equation for the conditional probability den- 
sity p(x, 1 1 Xq, s) in a time-dependent quadratic potential, 
U{x,t) = X[x - x min {t)] 2 /2 with x min (t) = (fi + f(t))/X, 



reading 

d t p = L(t)p = d x (U'(x,t)p) +Dd 2 xP , (2) 

with the absorbing boundary and initial conditions 

p(a, 1 1 xq , s) = for all t, s, and xq (3) 
p(x, s | x Q ,s) = 8(x - x ). (4) 

After firing the process immediately restarts at the in- 
stantaneous minimum of the potential. 

The set of eqs. defines our starting point for ob- 

taining the firing statistics of this driven neuron model. 
This is a rather intricate problem because the presence of 
non-stationarity and multiple time-scales for driving and 
relaxation, in combination with the absorbing boundary 
condition prohibits an analytical exact solution Our 
main objective is, nevertheless, to develop a most accu- 
rate analytical approximation that supersedes all prior 
attempts known to us. Those attempts, in fact, all in- 
volve the use of either of the following limiting approxi- 
mation schemes such as the limit of linear response theory 
(i.e. a weak stimulus f(t)) @, the limit of asymptotically 
weak noise 0,0 or the use of the method of images which 
appears to present an uncontrollable approximation for 
the case with A ^ 0, |]| . A most appealing numeri- 
cal approach is based on an exact integral equation for 
the first-passage time density of time-dependent Gauss- 
Markov processes with an absorbing boundary . Our 
scheme detailed below yields novel analytic and tractable 
expressions beyond the linear response and weak noise 
limit; these are limited solely by the use of a discrete, 
Markovian stochastic dynamics for the population of the 
attracting domain and slowly varying (in comparison to 
intra-well relaxation time-scale) stimuli f(t). As demon- 
strated below, this novel scheme indeed provides analyt- 
ical formulae that compare very favorably with precise 
numerical results of the full dynamics in eqs. Q] |21E}. 
Different from other approaches, we obtain the distribu- 
tion not only of the first-passage time but also of the 
residence time, which is the more interesting variable, 
concerning neurons. 
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To start, we approximate the solution to eqs. J2}0} 
in the regime where the statistics of times at which the 
threshold is reached can be characterized by a time- 
dependent firing rate rt(t) |^.Hl|. 

This rate then follows from a time-scale separation in 
the full FP dynamics J3J) with boundary condition J3J). 
After a few times A -1 of the fast relaxation the probabil- 
ity density p[x, t) assumes a slowly varying pattern that 
decays with the rate n(t). As in the time-independent 
case, this slowly varying part of p{x,t) can be expressed 
by a product of the (normalized) instantaneous station- 
ary solution po(x,t) oc exp{— U{x, t)/D} to the FP 
equation, satisfying L(t)po(x,t) — 0, and a form func- 
tion £(x,t). Our ansatz thus reads 



(t-»)>A -1 / t 

p(x,t | xq, s) ~ ((x, t)po(x, t) expl — Jk(s') ds' 



(5) 



The initial time s enters only through the exponential fac- 
tor. The dependence of the conditional probability on the 
initial value xq decays exponentially on the timescale A -1 
and therefore can be neglected for long times. 

Deep inside the attracting well, i.e. for x <C a there 
is no sensible difference between p{x,t) and po(x,t), and 
consequently £(x, t) approaches one. At the absorbing 
boundary, however, p(a,t) and consequently Ci a ^) both 
must vanish. The quantitative form of ((x, t) in the 
crossover region follows from 



L+(t)C(x,t)=0, 



(6) 



where the potential in the backward operator L + {t) = 
— U'{x, t)d x + Dd% can be linearized about the threshold 
x = a. Eq. (JSJ) then yields for the form function the result 



C(x, t) = 1 — exp 



{(. 



U'(a,t) 



(J) 



The rate n(t) is determined by multiplying the FP equa- 
tion J2J) in the long-time limit (J5J by the form func- 
tion £(x, t) and integrating over x from — oo to the thresh- 
old voltage a. In doing so, we account for prominent finite 
barrier corrections, yielding 

, f . _ J-oo dx CO. t)L(t)C(x, t)p (x, t) ^ 
f! oo dx( 2 (x,t)p (x,t) 

Upon insertion of eq. J7J) for the form function one can 
exactly perform the integrations and obtains for the rate 



«(t) = A 



AU(t) l-erf(y/AU(t)/D) 
D 1 -exp(-At/(t)/D) 



(9) 



where AU(t) denotes the instantaneous potential height 
at the threshold as seen from the minimum, and erf (z) is 
the error function. For very small D an expansion of the 



error function leads to the well-known weak noise result 
for the time-dependent rate 0,0] , i-e- 



K wn (t) = A y/AU(t)/(nD) exp(-AU(t)/D) . (10) 

Firing time distributions. — With the expression for the 
exit rate n(t) in @ we can calculate the properties of 
interest, namely the densities for the first-passage time 
and the residence time °f the attracting " integrating" 
state that covers the domain — oo < x(t) < a. 

The first-passage time distribution is given by the neg- 
ative rate of change of probability finding the process at 
time t in the "integrating" state, i.e. 



(t\s) = -d t / p(x, 1 1 x Q ,s) dx 



(11) 



= n(t) cxp(-j'K(s')dA , (12) 

s 

Here, the integral over the spatial part of p(x, t) yields 
unity and the exponential factor in J5| is obtained. It 
gives the probability for the process to stay in the "inte- 
grating" state from time s until t without interruption. 

The distribution of the residence times h(r), also 
termed the interspike interval density, follows as the aver- 
age of the first-passage time density over the density of re- 
setting times s, which coincides with the firing rate k(s). 
It thus reads: 



h(r) = lim 



)n(s)ds 



J_ T n(s)ds 



(13) 



Equations l|12|l and (|13|l . together with the expression 
for the rate © constitute the main results of this work. 
Their quantitative validity for an extended parameter 
regime will be checked next. 

Numerical comparison. — We have employed three 
methods for the numerical analysis. A first one is based 
on the Langevin equation where the position x(t) is 
updated sequentially. For the second we have solved the 
FP equation numerically, using a Chebychev colloca- 
tion method to reduce the problem to a coupled system 
of ordinary differential equations, see also hj. The third 
method solves an integral equation for the first-passage 
time probability density and is described in 0] . All three 
methods have provided practically identical results. 

In the prior stimulating work jj| it has been left open 
in what relevant parameter regime the employed approx- 
imation possesses validity p|. Here, using a periodic 
modulation f(t) = Acos(u>t) we like to determine a min- 
imal set of relevant parameters that can be taken for 
comparison with physiological measurements. The most 
general Langevin equation Q considered has seven con- 
stant parameters, including the threshold a and the reset- 
position xo. Three of them, (A, /i, a), can be chosen to 
be (1, 0, 1) by transforming to dimensionless coordinates, 
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using the time-unit A -1 , the space- unit (a — jit/A), an d 
the coordinate-origin at /i/A. The resulting parameters 
thus read (the bars indicate dimcnsionlcss coordinates) 

A = A/(Xa -fj,),w = lu/X, D = D/(X{a - fiX) 2 ) . (14) 

For our purposes it is advantageous to use instead the 
equivalent set 



U+/D, U-/D, andw, 



(15) 
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FIG. 1: First-passage time density for parameters U+/D — 8, 
U-/D = 5, u = 0.05. The jagged line depicts the his- 
togram obtained from iterations of the Langevin equation 
in 11611 . Note that fluctuations in the histogram depend on 
the total number of events and the width of the histogram 
bins. These fluctuations stay completely within their ex- 
pected range which is indicated as the gray shaded area in 
the inset. The height of this area is twice the expected stan- 
dard deviation of the histogram levels. The solid line shows 
the analytic first-passage time density from eq. 1121 , with the 
rate used in 10. Both lines are in excellent agreement. 
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FIG. 2: Residence time density vs. time for the same param- 
eters as in Fig. The jagged line shows the histogram ob- 
tained from iterations of 1161 . Again, the numerics practically 
coincides within the line-width with the analytic expression 
in I13H evaluated with the rate in © (solid line). 




where U+ and U— are the maximum and the minimum 
of U(l, t) during a whole period of modulation. These so 
chosen parameters have the benefit that they provide an 
on-hand estimate for the validity of our approximations 
and can be evaluated directly from (|TJ- 

The equation we have used in our simulations thus 
reads (omitting the overbars) 



x(t) = -x{t) + A cos(uit) 



(16) 



with the threshold located at x = 1. For obtaining the res- 
idence time, x has been reset into the minimum of U(x, t) 
immediately after firing. The Figs. ^ and [21 depict the 
probability densities of the first-passage and the resi- 
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FIG. 3: Testing extreme limits. First-passage time density 
for an extremely small lower barrier U-/D — 3. The remain- 
ing parameters are as in Fig. Q Langevin simulation results 
(jagged), analytical result in I12H with eq. © (solid), like- 
wise, with eq. 1101 (dashed). The inset compares the rate 
n(t) obtained from simulations of 1161 with the analytic re- 
sults from eq. ||UJ and eq. HOI , respectively. 
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FIG. 4: Testing extreme limits. Probability density of the 
first-passage time for a very fast driving, u — 0.5. The 
other parameters are as in Fig. Q The analytic approxi- 
mation (eqs. I12H and ©, dashed line) still depicts maxima 
that are approximately located at 1,2,... while the numeri- 
cal results are shifted to later times. The jagged line presents 
the Langevin- iterations from 116H . The solid curve presents 
within its linewidth both the numerical solution to the FP 
equation and to the integral equation from |Tcj|. 
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dence times, respectively. The residence time distribu- 
tion exhibits a less pronounced modulation in compari- 
son with the first-passage time distribution. Both analyt- 
ical expressions in (|12f> and (|13f) compare very favorably 
with the numerical results obtained by iteration of the 
Langevin equation Ijltjll ■ The remaining deviations in the 
two figures are of purely statistical nature (see the inset 
in Fig. ^) and can be diminished further by increasing 
the number of events in the simulations. 

In order to further test the range of validity of our 
novel approximation scheme we - on purpose - have cho- 
sen extreme values for the lower barrier height U-/D 
and angular driving frequency w, respectively, see Figs.|3] 
and Here, deviations from the numerical results are 
not the result of statistics but are systematic. For the low 
potential barrier in Fig. the time-scales in the process 
are not separated sufficiently. The fast intra-well fluctu- 
ations begin to influence the behavior of the modulated 
firing dynamics. Moreover, the difference between the 
moderate noise result for the time-dependent rate n(t) 
in (jjJJ and its weak noise approximation in (|10|l becomes 
visibly increased, as expected. Figure ^depicts the other 
extreme situation with a modulation time-scale that is 
not slow enough. Because the system cannot follow the 
driving instantaneously we find a shift in the maxima 
of the first-passage time density. This shift is not repro- 
duced by our approximation in and (fT2|> : nevertheless, 
our scheme yields amazingly good results even within this 
extreme parameter regime. The results based on the nu- 
merical evaluation of the FP equation J2J) and the integral 
equation 0] virtually collapse into one curve and per- 
fectly coincide with the Langevin simulations. The same 
result was obtained for the other parameter values. 

Conclusions. — The precise theoretical modeling of the 
neuronal spiking activity under external time-dependent 
driving presents a challenge of considerable importance 
in neurophysiology and physics. Due to the presence of 
non-stationarity, absorbing boundary conditions of the 
underlying first-passage problem and differing time-scales 
the task of obtaining reliable analytical estimates for 
the firing statistics is anything but trivial. By refer- 
ence to a discrete Markovian dynamics for the corre- 
sponding full space-continuous stochastic dynamics we 
succeeded in obtaining analytical approximations for the 
time-dependent first-passage time and the residence time 
statistics that are valid beyond the restraining limits of 
linear response and asymptotically weak noise. We have 
tested our findings for the case of a periodically driven 
LIF model. The obtained agreement with precise nu- 
merical simulations of either the Langevin type in 
or, equivalently, of the FP type in © turns out to be very 
good. Our method is not restricted to an oscillatory forc- 
ing but applies as well to arbitrary drive functions f(t) 
such as an exponentially decaying drive (e.g. simulat- 
ing a decaying threshold). Our scheme even yields good 
results in extreme parameter regimes where agreement 



cannot be expected a priori. 

Our method, primarily aimed at describing first- 
passage time and residence time probabilities of driven 
dynamical systems, is also readily extended to more re- 
alistic neuron models such as e. g. t he two-dimensional 
driven FitzHugh-Nagumo model 14] for neuronal spik- 
ing activity, whose multiple attractors may be consid- 
ered as discrete states. Likewise, the scheme can also 
be employed to study yet other time-dependent switch- 
ing dynamics and synchronization phenomena such as 
the paradigm of Stochastic Resonance la] and discrete or 
continuous Brownian motor transport |l5l ] . 
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